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Abstract. We outline the eccentricity evolution of sub-parsec massive black hole binaries 
(MBHBs) forming in galaxy mergers. In both stellar and gaseous environments, MBHBs are 
expected to grow large orbital eccentricities before they enter the gravitational wave (GW) 
observational domain. We re-visit the predicted eccentricities detectable by space based laser 
interferometers (as the proposed ELISA/NGO) for both environments. Close to coalescence, 
many MBHBs will still maintain detectable eccentricities, spanning a broad range from < 1CP 5 
up to <0.5. Stellar and gas driven dynamics lead to distinct distributions, with the latter 
favoring larger eccentricities. At larger binary separations, when emitted GWs will be observed 
by pulsar timing arrays (PTAs), the expected eccentricities are usually quite large, in the range 
0.01 — 0.7, which poses an important issue for signal modelling and detection algorithms. In this 
window, large eccentricities also have implications on proposed electromagnetic counterparts to 
the GW signal, which we briefly review. 



1. Introduction 

Multimessenger astronomy with massive black hole (MBH) binaries (MBHBs) has been discussed 
frequently in recent literature (e.g. [3-S]) as a unique way of fully understanding the mutual 
interaction of these systems with their environment. The general picture of MBHBs forming 
after a galaxy merger has been given in [§], but observational evidence for sub-parsec MBHBs 
is still lacking (see fiol ] for a review). Numerical studies down to parsec scales [ll- 14 ] have 
shown that the two MBHs continue to lose orbital energy and angular momentum under the 
large-scale action of gas/star-dynamical friction, and end up forming a Keplerian binary. At this 
point, the gaseous and stellar environment furthers the inspiral by either three-body scattering 
of individual stars and/or the interaction of the binary with a circumbinary gaseous disc (e.g. 
[HI Ell)- At this late evolutionary stage, excitation of e was already noticed and studied in both 
environments (see, e.g., @, [T3]). 

A non zero source eccentricity has important implications for gravitational wave (GW) obser- 
vations: on one hand, eccentricity has an impact on source detection and parameter estimation; 
on the other hand, measured eccentricities carry valuable information about the dynamical evo- 
lution of the detected systems. In gas rich environments, the surviving detectable eccentricity in 
the LISA band (10 -4 — 10 _1 Hz) was addressed in 0], stating that ens A ~ 10~ 3 — 10 _1 depending 
on the MBHB parameters. However, in view of the recent re-design of the LISA mission, esti- 



mates have to be re-adapted accordingly. The New Gravitational Observatory (ELISA/NGOj]] 
will see a shorter portion of the binary inspiral, with inevitable consequences for parameter 
estimation, and for residual eccentricity detectability. Pulsar timing arrays (PTAs) [l8|] are 
expected to detect the first GWs in the low-frequency band by the end of this decade. Given 
their frequency window (10~ 9 — 10~ 7 Hz), they will be sensitive to very massive (> 1O 8 M0) 
systems only , whose chirping is still marginal. This implies at the same time that observed 
MBHBs might still be coupled to their environment. We assess here the expected eccentricity 
distribution for these heavy sources in the PTA band. We then focus on the possibility of having 
observable PTA systems that are to some extent still coupled to either their stellar or gaseous 
environment, and assess the plausibility of detecting possible EM counterparts with current and 
up-coming missions, such as eROSITA \v^, MAXI Q and ATHENA!. 

The paper is structured as follows: First we describe the mechanisms of eccentricity growth 
in stellar environments in Sec. 12.11 and gaseous discs in Sec. 12.21 We then calculate the residual 
eccentricity distributions in the ELISA/NGO band in Sec. 13.11 and comment on the expected 
eccentricity population in the PTA band in Sec. 13.21 We review the implications this has on 
multimessenger signals in Sec. 13.21 and give our conclusions in Sec. 01 

2. Mechanisms of eccentricity growth: 

In the following we discuss the main physical mechanisms driving the MBHB eccentricity evo- 
lution in the dense astrophysical environment of a post-merger galaxy. We consider a MBHB 
with total mass M = Mi + M2 (-M2 < M\ are the individual MBH masses), semimajor axis o 
and eccentricity e. Throughout the manuscript, quantities without subscripts (such as M, a, e, 
E, L) always refer to the MBHB, whereas subscripts '1' and '2' refer to the individual MBHs. 

2.1. Stellar environment 

In dense stellar environments, the MBHB evolution is basically driven by three body interactions 



with individual stars [2114231]. To get a sense of the general picture, it is useful to consider the 
interactions between the stars and the binary one by one, as uncorrelated events, even though 
secular effects may also play a role (see, e.g., (23|). 

At the moment of pairing at sub-pc separations, the MBHB is expected to be surrounded 
by a dense stellar distribution of stars. If a is the MBHB semimajor axis, three body interac- 
tions with stars approaching the binary center of mass within a result in superelastic scattering, 
usually followed by the ejection of the stars. These stars carry away energy and angular mo- 
mentum, causing the binary to shrink. The efficiency of this mechanism depends on several 
environmental factors, such as the details of the stellar distribution and the effectiveness of 
relaxation mechanisms. We are here interested in the eccentricity evolution. 

Assessing the eccentricity evolution is not straightforward, because it depends on a com- 
bination of energy (E) and angular momentum (L) exchange during the scattering: Ae oc 
CiAE + CjAL, where the two c coefficients depend on the property of the MBHB. Let us con- 
sider, for simplicity, a bound star orbiting the MBHB at a given semimajor axis a*. The orbital 
energy that can be extracted in the scattering is oc a" 1 , whereas the angular momentum of the 
orbit is proportional to yfal. It follows that Ae oc — Qa^ 1 + Cjyfal (25[. If a* is small enough, 
then the first term dominates and the binary circularizes; vice versa, if a* is large enough, then 
the second term dominates, and the binary eccentricity grows. This behavior has to be expected, 
since stars with large a* have large angular momentum but just little energy, and the opposite is 
true for stars with small a*. It turns out that the transition point is at a* ~ a, i.e., stars within 
the binary orbit on average promote circularization, whereas stars outside the binary orbit make 



1 https : //lisa-light . aei .mpg. de/bin/view/ 

2 http : //www.mpe .mpg. de/athena/home .php?lang=en 



the binary grow eccentric. Heuristically this can be interpreted as follows. Let us consider a 
MBHB with M2 <C M\ and a given non-zero eccentricity. Such a binary spends most of its 
period near its apocenter, so in the case a* > a the probability of a close star-binary encounter 
(and subsequent star ejection) is maximal there. In the unequal mass case, the star is ejected 
following a series of close encounters with Mi. In such interactions the star gets a small 'kick' 
and the velocity of M2 decreases. As this velocity is perpendicular to the MBH position vector 
close to the apocenter, the secondary is forced on a more radial orbit. On the contrary, a star 
with a* < a "feels" the secondary hole when this approaches the primary at a distance < a. At 
that point, the interaction with the star is unlikely to occur close to the pericenter of the MBHB 
orbit (because the time spent by the binary at the pericenter is very small), and typically also 
extracts a large radial component from the velocity of the secondary black hole, hence causing 
circularization. Considering the collective effect of the stellar distribution, the evolution of the 
orbital eccentricity is therefore governed by the relative weight of stars inside and outside a. The 
steeper is the stellar cusp, the less effective is the eccentricity growth. As a general trend, quasi 
circular-equal mass MBHBs experience just a mild eccentricity growth, while systems which are 
already eccentric at the moment of pairing, or with mass ratio significantly lower than 1. can 



evolve up to e > 0.9 [17]. Such a trend has been observed in several numerical studies [111. 12614281] . 



The above argument is strictly valid for isotropic stellar distributions. The evolution of 



the binary eccentricity can be extremely different for non isotropic systems. For example, [291 ] 
demonstrated that in a rotating stellar system, the eccentricity evolution of unequal MBHBs is 
dramatically affected by the level of co/counter rotation of the stellar distribution with respect 
to the binary, with corotating distributions promoting circularization rather than eccentricity 
growth. This is because, as a simple consequence of angular momentum conservation during the 
ejection process, stars counterrotating with the binary tend to extract a lot of angular momentum 
from the MBHB, causing the eccentricity growth, whereas corotating stars do not. Nonetheless, 



most of the simulations involving rotating bulges [30j, |31[, or merging systems [32j, |33J find 
quite eccentric binaries at the moment of pairing (ranging from 0.4 to 0.8), and the subsequent 
evolution leads to a general eccentricity growth, in good agreement of what predicted for an 
isotropic stellar distribution. This may be because the binary evolution is mostly driven by loss 
cone refilling of unbound stars on almost radial orbits, with negligible initial angular momentum. 



On the other hand, [3J] find milder eccentricity evolutions, along the lines of what is expected 
for a rotating stellar remnant. The eccentricity evolution of MBHBs in stellar environments 
is still a largely open issue, but a significant eccentricity seem to be a general outcome of the 
models. 

2.2. Gaseous environments 

The study of the secular evolution of a MBHB immersed in a gaseous disc is essentially similar 
to type-II planet migration [3^]. However, when self-gravity of the disc and the accretion of 
matter onto the MBHs is included, the clear predictions of resonances and decay times coming 
from perturbation theory tend to be smeared out. In this paragraph, we will sketch the two 
effects that lead to the excitation and the saturation of the eccentricity evolution: 

• excitation of eccentricity via asymmetric torques from the circumbinary disc 

• saturation of eccentricity growth due to accumulated gas around the black holes 

When dynamical friction becomes inefficient, the subsequent MBHB evolution relies on transfer 
and redistribution of energy and angular momentum from the binary to the so-called circumbi- 
nary disc. This usually occurs at sub-parsec binary orbital separations, and goes hand in hand 
with the appearance of a so-called gap. A gap is a hollow, low-density region which is found to 
have a radial extent 5 of about twice the semi-major axis a of the binary (e.g. [g, 3^, 3?|]) which 



is attributed to the excitation of outer Lindblad resonances (OLR) with the most prominent 



one at 8 m =\ = 2 2//3 a ~ 1.59 a H. However, gas will still continue to leak from the inner ed ge o f 
the disc to the black holes and form mini-discs inside the respective Roche-lobes @, 0, [H,l3§] . 
As first shown by 0] and confirmed by [13], a binary inside a gap will not remain circular and 
neither will the inner edge of the disc (see also [36|, 0, |4l(] ) . 

As showed in [6] for a MBHB of total mass M and mass-ratio q = M2/-M1 = 1/3, the eccentricity 
of the binary will continue to grow up to about e ~ 0.6 while the eccentricity of the disc stays 
small, edisc < 0.15, even for e > 0.7. The reason for this growth of e can be understood by com- 
paring the angular velocity of the black holes to that of the fluid elements in the circumbinary 
disc. The movement of the MBHB inside the disc will induce overdensities close to the inner 
rim of the circumbinary disc, moving at lower angular speeds than the MBHs. If the MBHs 
are at apocenter, where they are slowest and thus remain the longest, they will experience a 
deceleration (negative torque) thus decreasing their associated angular momenta L\p and in- 
creasing e. Since the disc is at the same time expanding radially to accommodate the larger 
maximum MBHB separation, the theoretical limit to this growth depends on the evolution of 
8(e, t). Equating the average angular velocity of the inner rim of the disc to that of the binary 
at apocenter: (wfl u id) = oj &po yields the implicit relation for a limiting e: 

{ l±^ = 5 a (e,tf, 5 a = 5/a (1) 
(1 - e) 

which shifts the question to determining the function 5(e,t). Assuming S(e,t) = const = 2, 
0| found a limiting eccentricity of ~ 0.6. However we expect 5(e, t) to generally be a growing 
function of e related to the strongest OLRs. A detailed discussion of the OLRs and their 
evolution will soon be found in j42|. In general, if 5 grows with e, Lfl u id < -^2 is always true, and 
e continues to grow unless some other effect comes into play. When the binary is very eccentric, 
the Roche-lobes and the mini-discs attached to each individual MBH overlap at pericenter. 
Whenever this starts happening the secondary MBH will have to pass through the mini-disc of 
the primary and it will be slowed down at pericenter, decreasing e. This effect is thus strongly 
dependent on the thermodynamics of the gaseous disc and the mass contained within the mini- 
disc. 

In the remainder of this section, we will elucidate the growth of e with numerical examples 
of Newtonian self-gravitating discs around a MBHB with q = 1/3, using a setup described in 
the details of which will soon be found in [42 ]. The two examples use the same initial data 
and code, differing only in the thermodynamics treatment inside the gap. The dry run allows 
for an adiabatic evolution of the internal energy plus a /3 cooling inside the gap which effectively 
suppresses gas-inflows, whereas the wet run furthers the inflows by treating the gas isothermally 
when crossing the radius of r; so = 1.73 a, thus allowing us to differentiate the effects of the 
gas in the gap. When plotting the torque density T exerted by the disc onto the MBHs in 
figure [H we show the individual torques in the ^-direction onto each MBH in the middle and 
right panels and the sum T z = M\ T Zl + M2 T Z2 in the left panel. Since the binary is rotating 
counterclockwise, the negative portion of the torques always stems from the upstream region 
(yellow/red). Comparing the top and the bottom set of panels, we observe that altering the 
prescription for the thermodynamics inside the gap alters the torque density and the dynamical 
evolution of the entire disc. It is clear that the torque from the mini-discs are larger for the wet 
run (bottom panels), since the mass inside the mini-discs is higher and the gas is not allowed 
to cool (i.e. its internal energy remains constant) after crossing the threshold radius rj so . The 
torque density onto M\ (middle panel) is more or less symmetric in the azimuthal direction, 
whereas M2 clearly experiences larger torques from the part of the disc that is close to it. What 



3 The excitation of the m = 1 OLR is suppressed in equal mass binaries, nevertheless similar gap-sizes are also 
found for exactly equal mass MBHBs 36] 
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Figure 1. Torque density in the z-direction for two initially circular MBHB runs reaching 
e ~ 0.1 after 90 orbits. Top (bottom) panels: dry (wet) thermodynamic prescription for the 
cavity. In each line, the three panels are, from left to right: the total torque density on the 
binary T z = Mi T Zl + T Z2 , torque onto the primary MBH T Zl , torque onto the secondary 



MBH T Z2 . Figure made using SPLASH 43] 



is moreover notable, is the fact that in the left panel, where the total torque is shown to have 
a quadrupolar structure, the contributions from the mini-discs do not have exactly the same 
orientation as the contributions from the disc. However, the growth of e is very similar for 
wet and dry runs, which leads us to conclude that the excitation of e stems mostly from the 
circumbinary disc and is robust under different numerical prescriptions. The precise value of 
eiimit; however, depends on the density and angular momentum of the mini-discs, but is expected 
to be eiimit £ [0.6,0.8]. Further studies of the precise dependence on the thermodynamics are 
ongoing 421 ] . 



3. Implications for ELISA/NGO and PTA sources 

Low-frequency GWs are expected to be rich in astrophysical information if the source parameter 
analysis is sufficiently accurate to differentiate between various proposed astrophysical scenarios. 
Great hopes lie in a coincident detection of EM counterparts to complement the GW information 
with an independent measurement of the source redshift and a connection with the host galaxy, 
opening new scenarios for precision cosmology (see, e.g, [l|). The proposed ELISA/NGO instru- 
ment is designed for 10 3 — 10 7 M MBHBs chirping and merging in the 10" 4 - 10" 1 Hz band, 
whereas PTAs will be sensitive to GW frequencies in the range 10~ 9 — 10 _6 Hz with successful 
detections depending on the sensitivity of the radio instruments, the number of available stable 
millisecond pulsars and the observation time. 

To guide the development of effective detection pipelines, it is reasonable to review the as- 



sumptions of GW sources emitting purely sinusoidal, monochromatic signals. To estimate the 
impact of the eccentricity growth via the mechanisms described above, we need to understand at 
which binary separation they are effective. We therefore compare the star/gas driven shrinking 
timescales to the GW energy loss timescale. The GW timescale is given by: 

t GW = 7.84 x 10 7 yrM 8(?s - 1 ^F(e)- 1 , F(e) = (1 - e 2 )" 7 / 2 (l + ^e 2 + ge 4 ) (2) 

where q s = Aq/(l-\-q 2 ) is the symmetric binary mass ratio, Mg = M/10 8 M and is the binary 
separation in units of 10 3 Rs (Rs = 2GM/c 2 ). 

In gas rich environments, assuming a /3-disc model 4Jj for the circumbinary disc, the migra- 
tion timescale is 0] 



/ ■ \ -5/8 

t m = 2.09 x 10 6 yra-f {^f j M 8 3 V^(e) 7 / 8 af. (3) 

Here we defined the viscosity parameter ao.3 = a/0.3, the accretion efficiency eo.i = e/0.1 and 
the mass accretion rate rho.3 = m/0.3, where m = M/MEdd is the accretion rate normalized to 
the Eddington rate. 

In a star dominated system, assuming a full loss cone [ljj the MBHB hardening timescale is 
instead given by 

t h = 2.89 x 10 6 yrawoM^p^a^H^ 1 , (4) 

where we introduced the stellar velocity dispersion aioo = cr/100 Km s -1 , the stellar density 
at the binary influence radius (see [l7l ]. for details) ps = p/10 5 M pc -3 and the dimensionless 
hardening rate = H/15. Note that the combined M g _1 and 1 scaling implies that is 
independent on the MBHB mass. 

By equating t m and to i gw , one gets the characteristic a at which GW emission takes over 
in the binary evolution. It turns out that binaries are still coupled to their environment down to 
a separation of few hundred Schwarzschild radii, where the most massive systems fall in the PTA 
frequency band. This has two direct consequences: (i) PTA sources are generally eccentric, with 
an eccentricity distribution directly determined by the shrinking mechanism, and (ii) significant 
residual eccentricity might be present close to coalescence, in the ELISA/NGO band. 

In the following we apply selfconsistent models for the evolution of the MBHBs in gaseous 
and stellar environments to the expected population detectable by LISA/NGO and PTAs. In 
gas rich systems, we assume MBHBs to maintain an eccentricity of eii m it = 0.6 as long as they 
are coupled to their environment and then to circularize under the action of GW emission after 
decoupling. In star dominated systems, we use the self-contained model for the evolution of the 



MBHB under the effect of star ejection and GW radiation as outlined in 17] 



3.1. Eccentricity population of ELISA/NGO sources 

In view of the new sensitivity curve of the ELISA/NGO instrument, we re-visit here the ques- 
tion of the residual detectable eccentricity of a population of MBHBs evolving either in stellar 
or gaseous environments. We take the population of merging MBHBs predicted by a standard 
MBH cosmic evolution model where MBHs grow by merger and accretion starting from a pop- 
ulation of light seed MBHs, thought to be remnants of the PopIII stars [451 ] . To be specific, 
we took the model "LE" from [46J. We stress here that in this model, MBHBs were assumed 
to instantly coalesce after a dynamical friction timescale from the pairing of the two galaxies 
hosting the two MBHs: none of the dynamical evolution models is implemented selfconsistently 
in deriving the MBHB population. We therefore have a single MBHB population to which 
we apply different dynamical evolution models to determine the residual eccentricity. For each 




Figure 2. Eccentricity population of MBHBs detectable by ELISA/NGO and PTAs, expected 
in stellar and gaseous environments. Left panel: The solid histograms represent the efficient 
models whereas the dashed histograms are for the inefficient models. Right panel: solid his- 
tograms include all sources producing timing residuals above 3 ns, dashed histograms include 
all sources producing residual above 10 ns. 



mechanism (gas/star) we consider two scenarios (efficient /inefficient), to give an idea of the 
expected eccentricity range. The models are the following 

(i) gas-efficient: a = 0.3, m = 1. The migration timescale is maximized for this high values 
of the disc parameters, and the decoupling occurs in the very late stage of the MBHB 
evolution; 

(ii) gas-inefficient: a = 0.1, m = 0.1. Decoupling occurs at much larger (factor 3-to-5) separa- 
tions, and MBHBs have much more room to circularize; 

(iii) stars-efficient: p(ei) oc e. Initial eccentricities are taken according to a thermal probability 
density function [331 ] : 

(iv) stars-inefficient: e« = 0. Binaries are initially circular, a condition that minimizes the 
effectiveness of eccentricity growth. 

Results are shown in the left panel of figure El where we plot the residual eccentricity of all 
sources detected with SNR> 8, when they enter in the ELISA/NGO bancQ. In general, MBHBs 
evolving in gaseous environments are expected to retain a larger eccentricity, with a distribution 
tail possibly extending to e ~ 0.5. Residual eccentricities are mostly in the 10 -4 — 10 -2 range, 
and are generally lower for star driven systems. Note that the distribution predicted by the 
star-efficient model just overlaps to the one predicted by the gas-inefficient one. Therefore, 
apart from those two extreme cases, the distributions predicted by the two families of models 
are generally distinct, especially at the high e tail. A sufficient source detection statistics should 
allow discrimination between gas and star driven dynamics, providing valuable information of 
the MBHB environment. 



4 The band entrance is denned as the frequency where the strain of the source h s squared equals the one-sided 
noise spectral density of the detector S(f). 



3.2. Eccentricity population of PTA sources 

To investigate PTA sources, we instead use MBHB population models from (4?l |. Here we 
just want to highlight the fact that PTA sources might be very eccentric, and adequate signal 
processing algorithms need to be developed. We therefore use a single representative model for 
each mechanism. For the gas driven evolution we use a = 0.3, m = 0.3; whereas for the star 
driven model we use = 0.1 for all the binaries. In the right panel of figure [21 we plot the 
eccentricity distribution of all sources contributing to the signal at timing residual level above 
10ns or 3ns, which is expected to be in the reach of future PTA campaigns. The distributions 
are rather similar with more spread and higher extremal eccentricities for the stellar case. It 
is not clear, whether the two distributions are discernible; the figures, however, clearly suggest 
that a large portion of PTA signals would have fairly high eccentricities (e > 0.1). 



Multimessenger signals. Concentrating on gaseous disc environments, [7[ and [8] have, in a 
complementary way, outlined possible counterparts to resolvable PTA sources, which we will 
quickly review here. First, we will look at coupled systems, as done in 0], then conclude with 
a remark on how eccentricity impacts on decoupled systems as studied in [8j. 
If t v is the viscous time of the disc, then as long as tow > t u , the disc can dissipate inwards fast 
enough to follow the MBHB and during this phase the MBHB is expected to be highly eccentric. 
In [6], it has been shown that the strength of a periodicity in the light curve related to the MBHB 
orbital frequency depends on the eccentricity of the orbit and the angular momentum distribution 
of the inspiralling gas. First, it is important to note that for circular or mildly eccentric orbits 
the periodicity in the gas inflow rate is much weaker than for high (e >0.4) eccentric systems. 
The number and quality of observations needed to identify the orbital periodicity is challenging 
for upcoming all-sky surveys, even assuming highly modulated inflow rates for e ~ 0.6 systems. 
This is simply because of the high number of pointings per orbit needed (which will be the main 
limitation of eROSITA) and the required flux limit of the instrument (which limits ongoing 
surveys like MAXI). Detecting a PTA counterpart in a gaseous environment will thus be likely 
biased towards eccentric systems that accrete the instreaming gas on time scales comparable 
to the binary period. If, however, in the above scenario the time scale for the accretion of 
the mini-disc is long compared to the binary period, the instreaming gas periodicity might be 
diluted in the process of angular momentum redistribution within the disk, and might not be 
reflected in the accretion luminosity onto the two MBHs. In this case, stable, broad double Ka 
lines arising from the inner portion of the two mini-discs might be detectable by ATHENA 0], 
and could provide meaningful measurements of the two MBH spins. If, on the contrary, the gas 
is instantaneously accreted, there will be no detectable lines and the general contribution from 
of gap region to the emitted luminosity will likely be negligible; only if the spread in inflo wing 
angular momentum is very large, colliding caustics can form small inviscid accretion discs [48l. l49lj 
making the MBHs variable, luminous X-Ray sources. 

Looking at circular binaries, [s| have shown that gas can continue to leak towards the MBHB 
even after decoupling, finding surface densities as large as £ ~ 10 4 g/cm 2 at ~ 40 GM/c 2 for a 
total MBHB mass of 10 9 M©. While for eccentric systems the width of the gap is generally larger 
than for circular ones, their Green function approach is valid nonetheless, as long as the leaking 
gas is ripped out of the edge of the circumbinary disc on dominantly radial orbits, i.e. having a 
much lower angular momentum than the inner edge of the circumbinary disc. Their results (see, 
e.g. their figure 4) would not change much qualitatively, if they instead used a larger cavity to 
begin with and then adapted the inner boundary condition to a function 5 a (e(t)) rather than 
using 2a(t) . However, due to gas from the decoupled disc being required to diffuse over larger 
distances while the gas is also colder at the same time thus less diffusive, it is expected that the 
total amount of fossil gas is smaller the higher the eccentricity. 



4. Conclusions 

We have shown that the excitation of MBHB eccentricity in both star and gas dominated 
environments is a robust prediction of state of the art dynamical evolution models. In the 
gas dominated case, limiting eccentricity values eii m it > 0.6 are expected following excitation 
by a cold dense circumbinary disk, whereas stellar dynamics can easily produce eccentricities 
higher than 0.9. As a consequence, PTA sources, many of which are at the evolutionary stage 
in which the environment still plays a role in their dynamical evolution, are expected to be 
significantly eccentric. Many binaries producing residuals at a 10ns level (in the reach of future 
PTA campaigns) have e > 0.1. The impact of such high eccentricities on the PTA signal should 
be appropriately taken into account when developing dedicated data analysis algorithms. For 
ELISA/NGO sources, which have long since decoupled from their environment, we have shown 
the distribution of the residual eccentricity expected for both environments. As the latter 
are disjoint at the high-eccentricity tail, enough GW observation statistics could distinguish 
among the different scenarios, providing valuable information on the environment of coalescing 
MBHBs. As far as multimessenger signals are concerned, we did not treat electromagnetic 
counterparts to ELISA/NGO events and focused on PTA observations. In this latter case, the 
electromagnetic signal is generally bright, due to the high mass and relatively small luminosity 
distance of the sources, making certain distinct features detectable by upcoming observatories. 
We have shown that counterpart identification is likely to be bias towards eccentric MBHB in 
periodicity searches, whereas, depending on the angular momentum of the instreaming material 
feeding the mini-discs, peculiar double Ka lines might be identifiable by future hard X-ray deep 
spectroscopy. 
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